﻿void SearchLinkCellList(caja *bx)
{
	int c, c1;
	int iii,jjj;
	int shift[3];

	for(int i=0;i<bx->npart;i++)
	{
		bx->part[i].nvecinos=0;
	}

	//Recorremos todas las celdas
	for (int ix=0;ix<bx->ncell;ix++)
	{
		for (int iy=0;iy<bx->ncell;iy++)
		{
			for (int iz=0;iz<bx->ncell;iz++)
			{
				///Calculamos el indice de la celda (ix,iy,iz)
				c = int(ix*bx->ncell2+iy*bx->ncell+iz);
				///Recorremos las particulas de esta celda y las celdas vecinas
				for (int ixred=ix-1;ixred<=ix+1;ixred++)
				{
					/*shift: condiciones de contorno periodicas tenidas en cuenta para
					localizar las celdas vecinas*/
					if (ixred < 0)
					{
						shift[0] = int(ixred+bx->ncell);
					}
					else
					{
						if (ixred>=bx->ncell)
						{
							shift[0] = int(ixred-bx->ncell); 
						}
						else
						{
							shift[0] = int(ixred);
						}
					}
					for (int iyred=iy-1;iyred<=iy+1;iyred++)
					{
						if (iyred < 0)
						{
							shift[1] = int(iyred+bx->ncell);
						}
						else
						{
							if(iyred>=bx->ncell)
							{
								shift[1] = int(iyred-bx->ncell);
							}
							else
							{
								shift[1] = int(iyred);
							}
						}
						for (int izred=iz-1;izred<=iz+1;izred++)
						{
							if (izred < 0)
							{
								shift[2] = int(izred+bx->ncell);
							}
							else
							{
								if(izred>=bx->ncell)
								{
									shift[2] = int(izred-bx->ncell);
								}
								else
								{
									shift[2] = int(izred) ;
								}
							}
							///Calculamos el indice c1 de la celda vecina
							c1 = int(((shift[0]+bx->ncell)%bx->ncell)*bx->ncell2+((shift[1]+bx->ncell)%bx->ncell)*bx->ncell
								+((shift[2]+bx->ncell)%bx->ncell));
							///Buscamos la primera particula iii de la celda c
							iii = bx->head[c];

							while (iii != -1) ///iii=-1 ->No hay particulas en la celda
							{
								///Buscamos la primera particula jjj de la celda c1
								jjj =bx->head[c1];
								
								while (jjj != -1) ///jjj=-1->No hay particulas en la celda
								{
									//Evitamos el conteo doble de las interacciones
									if (iii < jjj) 
									{
										//La particula iii vecina de la part jjj
										bx->part[iii].lstVecinos[bx->part[iii].nvecinos]=jjj;
										bx->part[iii].nvecinos++;
									}
									//Pasamos a la siguiente particula de la celda c1
									jjj = bx->link[jjj];
								}
								//Pasamos a la siguiente particula de la celda c
								iii = bx->link[iii];
							}
						}
					}
				}
			}
		}
	}
}

double interaccionLinkCellList(caja *bx, bool calcStress)
{
	//Si calcStress = true, se calcula la parte correspondiente al potencial
	//del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
	int c, c1;
	int iii,jjj;
	int shift[3];
	int count=0;//Numero de pares interactuantes
	double dist,dist2,dist3,rm4,potdm, ePoten;
	double Co,Co2, angPart1,angPart2;
	double fx,fy,fz;
	double rc = bx->rcut;
	double rc2=rc*rc;
	double virial=0;
	//Tensor de esfuerzos
	double sxy, sxz, syz, szz;
	double rep;
	double dummy,f2;
	void stress(caja *bx);
	//Ponemos las fuerzas de las particulas de la caja a 0
	bx->initForces();

	//Inicializamos la energia potencial
	ePoten=0.;
	//Inicializamos el tensor de esfuerzos
	sxy=sxz=syz=szz=0.0;

	
	//Recorremos todas las celdas
	for (int ix=0;ix<bx->ncell;ix++)
	{
		for (int iy=0;iy<bx->ncell;iy++)
		{
			for (int iz=0;iz<bx->ncell;iz++)
			{
				///Calculamos el indice de la celda (ix,iy,iz)
				c = int(ix*bx->ncell2+iy*bx->ncell+iz);
				///Recorremos las particulas de esta celda y las celdas vecinas
				for (int ixred=ix-1;ixred<=ix+1;ixred++)
				{
					/*shift: condiciones de contorno periodicas tenidas en cuenta para
					localizar las celdas vecinas*/
					if (ixred < 0)
					{
						shift[0] = int(ixred+bx->ncell);
					}
					else
					{
						if (ixred>=bx->ncell)
						{
							shift[0] = int(ixred-bx->ncell); 
						}
						else
						{
							shift[0] = int(ixred);
						}
					}
					for (int iyred=iy-1;iyred<=iy+1;iyred++)
					{
						if (iyred < 0)
						{
							shift[1] = int(iyred+bx->ncell);
						}
						else
						{
							if(iyred>=bx->ncell)
							{
								shift[1] = int(iyred-bx->ncell);
							}
							else
							{
								shift[1] = int(iyred);
							}
						}
						for (int izred=iz-1;izred<=iz+1;izred++)
						{
							if (izred < 0)
							{
								shift[2] = int(izred+bx->ncell);
							}
							else
							{
								if(izred>=bx->ncell)
								{
									shift[2] = int(izred-bx->ncell);
								}
								else
								{
									shift[2] = int(izred) ;
								}
							}
							///Calculamos el indice c1 de la celda vecina
							c1 = int(((shift[0]+bx->ncell)%bx->ncell)*bx->ncell2+((shift[1]+bx->ncell)%bx->ncell)*bx->ncell
								+((shift[2]+bx->ncell)%bx->ncell));
							///Buscamos la primera particula iii de la celda c
							iii = bx->head[c];

							while (iii != -1) ///iii=-1 ->No hay particulas en la celda
							{
								///Buscamos la primera particula jjj de la celda c1
								jjj =bx->head[c1];
								
								while (jjj != -1) ///jjj=-1->No hay particulas en la celda
								{
									//Evitamos el conteo doble de las interacciones
									if (iii < jjj) 
									{
										///Calculamos la distancia con condiciones periodicas entre iii y jjj
										dist2= bx->distancia2(iii,jjj);
										//Si la distancia es menor que el cutoff calculamos la interaccion
										if (dist2 < rc2)
										{
											dist=sqrt(dist2);
											dist3=dist*dist*dist;
											rm4=1./(dist3*dist);
											count++;
											if (dist<4.0)
											{
												dummy=bx->u0*exp(-50.*(dist-2.));
												//Para las fuerzas
												rep=-0.1875*dummy;
												//Para el potencial
												f2=-0.00375*dummy;
											}
											else
											{
												rep=0.;
												f2=0.0;
											}

											///Calculamos la fuerza entre iii y jjj
											Co=bx->dr[2]/dist;
											Co2=Co*Co;																	
			
											//double rep = 0.1875*(bx->u0)*exp(-50.0*(dist-2.0));
											angPart1= (bx->u30)*(5.0*Co2-1.0)*rm4+rep;
											angPart2= (bx->u30)*(5.0*Co2-3.0)*rm4+rep;
											fx=angPart1*(bx->dr[0])/dist;
											fy=angPart1*(bx->dr[1])/dist;
											fz=angPart2*(bx->dr[2])/dist;

											if(calcStress)
											{
												sxy+=(bx->dr[0])*fy;
												sxz+=(bx->dr[0])*fz;
												syz+=(bx->dr[1])*fz;
												szz+=(bx->dr[2])*fz;
											}

											bx->part[iii].force[0]+=fx;
											bx->part[iii].force[1]+=fy;
											bx->part[iii].force[2]+=fz;

											bx->part[jjj].force[0]-=fx;
											bx->part[jjj].force[1]-=fy;
											bx->part[jjj].force[2]-=fz;

											potdm = (bx->u0)*(1.0-3*Co2)/(dist3);
											//Acualizamos la energia potencial
											ePoten+=2.*f2+potdm;
											//Actualizamos el virial
											//virial-=ccel3*(dr[0]*dr[0]-dr[1]*dr[1]-dr[2]*dr[2]);
											virial-=dist*(fx*bx->dr[0]+fx*bx->dr[1]+fx*bx->dr[2]);
										}
									}
									//Pasamos a la siguiente particula de la celda c1
									jjj = bx->link[jjj];
								}
								//Pasamos a la siguiente particula de la celda c
								iii = bx->link[iii];
							}
						}
					}
				}
			}
		}
	}

	if(calcStress)
	{
		//Actualizamos la parte del potencial del tensor de esfuerzos microscopico
		bx->spxy=sxy;
		bx->spxz=sxz;
		bx->spyz=syz;
		bx->spzz=szz;
	}
	//Introducimos la correccion del cutoff del potencial
	//stress(bx);
	bx->epoten=ePoten/bx->npart;
	bx->virial=virial;
	return ePoten;
}

double interaccionSimple(caja *bx, bool calcStress)
{
	//Si calcStress = true, se calcula la parte correspondiente al potencial
	//del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
	int c, c1;
	int iii,jjj;
	int shift[3];
	int count=0;//Numero de pares interactuantes
	double dist,dist2,dist3,rm4,potdm, ePoten;
	double Co,Co2, angPart1,angPart2;
	double fx,fy,fz;
	double rc = bx->rcut;
	double rc2=rc*rc;
	double virial=0;
	//Tensor de esfuerzos
	double sxy, sxz, syz, szz;
	double rep;
	double dummy,f2;
	void stress(caja *bx);
	//Ponemos las fuerzas de las particulas de la caja a 0
	bx->initForces();

	//Inicializamos la energia potencial
	ePoten=0.;
	//Inicializamos el tensor de esfuerzos
	sxy=sxz=syz=szz=0.0;

	
	//Recorremos todos los pares de particulas
	for(int iii=0;iii<bx->npart-1;iii++)
	{
		for(int jjj=iii+1;jjj<bx->npart;jjj++)
		{
			dist2= bx->distancia2(iii,jjj);
			dist=sqrt(dist2);
			dist3=dist*dist*dist;
			rm4=1./(dist3*dist);
			count++;
			if (dist<4.0)
			{
				//Si la distancia es menor que el cutoff calculamos la interaccion de esfera dura
				dummy=bx->u0*exp(-50.*(dist-2.));
				//Para las fuerzas
				rep=-0.1875*dummy;
				//Para el potencial
				f2=-0.00375*dummy;
			}
			else
			{
				rep=0.;
				f2=0.0;
			}

			///Calculamos la fuerza entre iii y jjj
			Co=bx->dr[2]/dist;
			Co2=Co*Co;																	

			//double rep = 0.1875*(bx->u0)*exp(-50.0*(dist-2.0));
			angPart1= (bx->u30)*(5.0*Co2-1.0)*rm4+rep;
			angPart2= (bx->u30)*(5.0*Co2-3.0)*rm4+rep;
			fx=angPart1*(bx->dr[0])/dist;
			fy=angPart1*(bx->dr[1])/dist;
			fz=angPart2*(bx->dr[2])/dist;

			if(calcStress)
			{
				sxy+=(bx->dr[0])*fy;
				sxz+=(bx->dr[0])*fz;
				syz+=(bx->dr[1])*fz;
				szz+=(bx->dr[2])*fz;
			}

			bx->part[iii].force[0]+=fx;
			bx->part[iii].force[1]+=fy;
			bx->part[iii].force[2]+=fz;

			bx->part[jjj].force[0]-=fx;
			bx->part[jjj].force[1]-=fy;
			bx->part[jjj].force[2]-=fz;

			potdm = (bx->u0)*(1.0-3*Co2)/(dist3);
			//Acualizamos la energia potencial
			ePoten+=f2+potdm;
			//Actualizamos el virial
			virial-=dist*(fx*bx->dr[0]+fx*bx->dr[1]+fx*bx->dr[2]);

		}
	}

	if(calcStress)
	{
		//Actualizamos la parte del potencial del tensor de esfuerzos microscopico
		bx->spxy=sxy;
		bx->spxz=sxz;
		bx->spyz=syz;
		bx->spzz=szz;
	}
	bx->epoten=ePoten/bx->npart;
	bx->virial=virial;
	return ePoten;
}

//double interaccionLinkCellList2(caja *bx, bool calcStress)
//{
//	//Si calcStress = true, se calcula la parte correspondiente al potencial
//	//del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
//	int c, c1;
//	int iii,jjj;
//	int shift[3];
//	int count=0;//Numero de pares interactuantes
//	void SearchLinkCellList(caja *bx);
//	double dist,rm4,potdm, ePoten;
//	double Co,Co2, angPart1,angPart2;
//	double fx,fy,fz;
//	double rc = bx->rcut;
//	double rc2=rc*rc;
//	double virial=0;
//	//Tensor de esfuerzos
//	double sxy, sxz, syz;
//	double rep;
//	double dummy,f,f2;
//	//Ponemos las fuerzas de las particulas de la caja a 0
//	bx->initForces();
//
//	//Inicializamos la energia potencial
//	ePoten=0.;
//	//Inicializamos el tensor de esfuerzos
//	sxy=sxz=syz=0.0;
//	
//
//	//
//	SearchLinkCellList(bx);
//
//	for(int iii=0;iii<bx->npart;iii++)
//	{
//		for(int j=0;j<bx->part[iii].nvecinos;j++)
//		{
//			jjj=bx->part[iii].lstVecinos[j];
//			///Calculamos la distancia con condiciones periodicas entre iii y jjj
//			dist= bx->distancia(iii,jjj);
//			//Si la distancia es menor que el cutoff calculamos la interaccion
//			if (dist < rc)
//			{
//				count++;
//				if (dist<4.0)
//				{
//					dummy=bx->u0*exp(-50.*(dist-2.));
//					//Para las fuerzas
//					f=-0.1875*dummy;
//					//Para el potencial
//					f2=-0.00375*dummy;
//				}
//				else
//				{
//					f=0.;
//					f2=0.0;
//				}
//
//				///Calculamos la fuerza entre iii y jjj
//				Co=bx->dr[2]/dist;
//				Co2=Co*Co;
//				rm4=pow(dist,-4.0);																			
//				rep=f;
//
//				//double rep = 0.1875*(bx->u0)*exp(-50.0*(dist-2.0));
//				angPart1= (bx->u30)*(5.0*Co2-1.0)*rm4+f;
//				angPart2= (bx->u30)*(5.0*Co2-3.0)*rm4+f;
//				fx=angPart1*(bx->dr[0])/dist;
//				fy=angPart1*(bx->dr[1])/dist;
//				fz=angPart2*(bx->dr[2])/dist;
//
//				if(calcStress)
//				{
//					sxy+=(bx->dr[0])*fy;
//					sxz+=(bx->dr[0])*fz;
//					syz+=(bx->dr[1])*fz;
//				}
//
//				bx->part[iii].force[0]+=fx;
//				bx->part[iii].force[1]+=fy;
//				bx->part[iii].force[2]+=fz;
//
//				bx->part[jjj].force[0]-=fx;
//				bx->part[jjj].force[1]-=fy;
//				bx->part[jjj].force[2]-=fz;
//
//				potdm = (bx->u0)*(1.0-3*Co2)/(dist*dist*dist);
//				//Acualizamos la energia potencial
//				ePoten+=2.*f2+potdm;
//				//Actualizamos el virial
//				//virial-=ccel3*(dr[0]*dr[0]-dr[1]*dr[1]-dr[2]*dr[2]);
//				virial-=dist*(fx*bx->dr[0]+fx*bx->dr[1]+fx*bx->dr[2]);
//			}
//		}
//	}
//
//	
//	if(calcStress)
//	{
//		//Actualizamos la parte del potencial del tensor de esfuerzos microscopico
//		bx->spxy=sxy;
//		bx->spxz=sxz;
//		bx->spyz=syz;
//	}
//	//Introducimos la correccion del cutoff del potencial
//	//ePoten-=LJcut*count;
//	bx->epoten=ePoten/bx->npart;
//	bx->virial=virial;
//	return ePoten;
//}

//void stress(caja *bx)
//{
//	double dr[3];
//	double sxy,sxz,syz;
//	sxy=sxz=syz=0.0;
//
//	for(int i=0;i<bx->npart-1;i++)
//	{
//		for(int j=i+1;j<bx->npart;j++)
//		{
//			for(int k=0;k<3;k++)
//			{
//				dr[k]=bx->part[j].pos[k]-bx->part[i].pos[k];
//			}
//			sxy+=dr[0]*(bx->part[i].force[1]);
//			sxz+=dr[0]*(bx->part[i].force[2]);
//			syz+=dr[1]*(bx->part[i].force[2]);
//
//		}
//	}
//	bx->spxy=sxy;
//	bx->spxz=sxz;
//	bx->spyz=syz;
//}

